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Modelling injection and feedback of Cosmic Rays in grid-based 
cosmological simulations: effects on cluster outskirts. 
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ABSTRACT 

We present a numerical scheme, implemented in the cosmological adaptive mesh refinement 
code ENZO, to model the injection of Cosmic Ray (CR) particles at shocks, their advection and 
their dynamical feedback on thermal baryonic gas. We give a description of the algorithms and 
show their tests against analytical and idealized one-dimensional problems. Our implemen- 
tation is able to track the injection of CR energy, the spatial advection of CR energy and its 
feedback on the thermal gas in run-time. This method is applied to study CR acceleration and 
evolution in cosmological volumes, with both fixed and variable mesh resolution. We compare 
the properties of galaxy clusters with and without CRs, for a sample of high-resolution clusters 
with different dynamical states. At variance with similar simulations based on Smoothed Par- 
ticles Hydrodynamics, we report that the inclusion of CR feedback in our method decreases 
the central gas density in clusters, thus reducing the X-ray and Sunyaev-Zeldovich effect from 
the clusters centre, while enhancing the gas density and its related observables near the virial 
radius. 
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1 INTRODUCTION 

Galaxy clusters form through the hierarchical assembly of dark 
matter (DM) and gas over cosmological times. This process hap- 
pens via super-sonic accretion, with efficient driving of turbulent 
motions, shock waves and magn etic field compress i on and am- 
plific a tion in the cosmic g as (e.g. iBvkov et al j|200ol ; iDolag et"ail 
l2008t iBriiggen et alj I20TII and references therein). This picture 



(e.g. iBellll 1971: iBlandford & Ostrikedl 19781: iDrurv & Voelklll98ll: 



Elliso n et alj Il995l: iKang & Jonesl 1 199(1 llvlalkov & O'C Drurv 



is based on observational evidences ( e.g. [Carilli & Tavloi l2002t 
iGovoniet al]|2002l ; iFerrari et alj|2008t for recent results ), and is 
theoretically supported by numerical simulations (e.g. Ryu et all 



1998llDolag et al.ll999l.llvliniati et al.l2 000, 200l i:lRvuet all 20031: 



Pfrommer et alj 120071; Hapichino & Niemevej [20081: iHoeft et all 
20081: IKang et alj 120071: Iskillman et alj l2008t lyazza et alj 1200*91 
Donnert et all 120091 : IVazza et alj I2009L l201ld : iBonafede et ail 
201 lb . However, the details of the interplay between various non- 
thermal components of the ICM are still quite uncertain. 

Shock waves driven by the accretion of matter constitute the 
most important channel for the thermalization of baryonic gas over 
cosmological time, and they may also be efficient sites of cos- 
mic ray particles (CR) acceleration. CRs are thought to be accel- 
erated in the pattern of MHD fluctuations across the shock tran- 
sition, in a process known as diffusive shock acceleration, DSA 
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I2001L IKang & Jonesll2007l : IKang et al.ll2009l ; ICaprioli et al]|201C 
The diffusion of accelerated particles across the shock produces 
a shock-precursor, which modifies the overall compression pro- 
duced by the entire shock structure in a non-linear way (e.g. 
IDrurv & Voelk1ll98ll : lBlasill2004l : lAmato & Blasill2005l> . The main 
factors which govern the overall efficiency of the injection of rel- 
ativistic particles are the Mach number of the shock, M, and the 
properties of the magnetic field in the shock transition. At quasi 
parallel shocks, M is the primary parameter which determines the 
acceleration efficiency of CRs, while the secondary parameter is 
the minimum momentum of injection of particles, p t h. required to 
start the acceleration. In the test-particle regime, where the acceler- 
ated particles do not have any dynamical feedback on the structure 
of the shock, the DSA predicts the universal energy spectrum of 
N(p)dp oc p~ q dp for the accelerated particles, where q is related 
to the Mach number only via q = 2(M 2 + 1) /(M 2 - 1). However, 
models of DSA predict that at strong shocks, M > 5, the fraction 
of accelerated CR particles becomes non-negligible compared to 
the thermal pool, £ > 10 -3 , and that a significant fraction of the 
shock energy is transferred to the CR particles, which in turn mod- 
ify the shock structure. The result is a larger shock compression and 
an amplification of the magnetic fields in the pre-shock region (e.g 
IKang et ai] [2009. for a modern review). The basic predictions of 
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DSA have been tested successfully with multi-band obse rvations of 
radio and 7 emission from remnants of su pernovae (e.g. lRevnoldsl 
l2008l:IVink et alj20ld : lEdmon et al.ll201 lh . 

Onc e accelerated, CR par ti cles can accumulate in galaxy 
clusters feerezinskv et al.l 1 1991 IVQlk & Atovan] [l999) possibly 
producing an important non-thermal component which could 
be directly sampled by future gamma ray observ ations (e.g. 
IPinzke & Pfrommej|201ol ; lBrunetti & Lazarianll201 lal) . Secondary 
particles injected in the ICM via proton- proton collisions may 
cause detectable synchrotron radiation (e.g Blasi & Colafrancesco 
ll999l:lDolag & EnBlinll2000h and, as well as primary particles, that 
are re- accelerated by MHD turbulence and may thus escape radio 
halos {jBrunetti et alj|2008l) . CR particles can mix with the thermal 
plas ma, and may have an important dynamical effect on the ICM 
(e.g. lRuszkowski & Ohll2"oTTI ; lBrunetti & Lazarianll201 lbh . Its role 
can also be that of shaping the evolution of X-ray cavities powered 
by AGN jets (e.g. iMafhews & Brighentil 120071 ; Iguo & Ohll20o3 ; 
Mathe ws &Guol201lh . and to provide an additional source of ther- 
malization in co ol core clusters, via Alfven wave dissipation (e.g. 
iGuo & Ohll2008h . 

The goal of this work is to provide a simplified but robust 
algorithm to model the run-time injection of CR particles at cos- 
mological shock waves, their advection and their dynamical ef- 
fect on the evolution of thermal gas in larg e-scale structures. This 
com plements the meth ods developed by Pfrommer et alj d2007l) 
and lEnfilin et alj ( 120070 , where the injection and feedback from 
CRs is simulated with Lagrangian cosmological simulations using 
smoothed-particle-hydrodynamics. To this end, we implemented 
new modules to follow CRs in the framew ork of the public 1.5 
version of the ENZO code l lBrvan et al.l 19951) . and we run a number 
of validation tests and large-scale production runs in cosmology. 

The paper is organized as follows: in Sect[2]we describe the al- 
gorithms developed to treat the injection, advection and dynamical 
feedback of CRs in ENZO; in Sect(3]we present one-dimensional 
tests for the code and its various options; in Sect|4] we present 
our 3-D results in cosmology, obtained with fixed mesh resolution 
(Sect l4.lt and adaptive mesh refinement (Sect l4.2t . A discussion of 
the results is given in Sect[5] while our conclusion are summarized 
in Sect|6] 



2 NUMERICAL METHODS 
2.1 ENZO 1.5 

All simulations in this work were performed using the ENZO 1.5 
code developed by the Laboratory for Computational Astrophysics 
at the University of California in San Diego (http://lca.ucsd.edu). 

ENZO is an adaptive mesh refinement (AMR) cosmologi- 
cal hy brid code highly optimized for high performance comput- 
ing (lBrvaneTaT]|l995l : lO'Shea et all |2004| ; iNorman et ai] |2007| ; 
ICollins et alJ2010h . It uses a particle-mesh N-body method (PM) to 
follow the dynamics of the collision -less Dark Matter (DM) com- 
ponent dHocknev & Ea stwood 1988), a nd an adap tive mesh method 
for ideal fluid-dynamics (Berger &Colelralll989h . 

The DM component is coupled to the baryonic matter (gas) 
via gravitational forces, calculated from the total mass distribu- 
tion (DM+gas) solving the Poisson equation with an FFT based 
approach. The gas component is described as a perfect fluid and 
its dynamics is calculated by solving conservation equations of 
mass, energy and momentum over a computational mesh, using an 
Eulerian solver based on the Piecewise Parabolic Method (PPM, 



IColella & Woodward|[T984h . This scheme is a higher-order exten - 
sion of Godunov's shock capturing method dGodunov et aljfl976l) . 
and it is at least second-order accurate in space (up to the fourth- 
order in 1-D, in the case of smooth flows and small time-steps) and 
second-order accurate in time. 

For a more detailed overview of the numerical algorithms in 
ENZO and of the updated strategies for handling the data-structures 
in the s imulations, we refer the readers to the excellent recent re- 
view of lCollins etai1d2010l) . 

2.2 Shock identification 

We follow the standard assumption of Diffusive Shock Accelera- 
tion, relating the injection efficienc y of Cosmic Rays at shocks to 
the Mach number of the shock (e.g. lKang & Jonesll2007L and refer- 
ences therein). 

In order to model the injection of CRs at shock waves, a reli- 
able on-the-fly method to measure the Mach number of shocks in 
the simulation is required. 

A few grid-based methods in the literature were tested in 
recent years that apply the Rankine-Hugoniot j umps conditions 
across shocked cells (e.g. iMiniati et alj I2OO0I: iRvu et alj 120031: 
ISkillman et all2008l ; lvazza et alj2009l ). 

In our algorithm, we detect and measure the Mach number of 
shock waves during run-time using an approach based on the dif- 
ferences of the gas pressure (or gas plus CR) between cells. This 
method is conceptually ve ry similar to the "Temperature Jump" 
method of |Rvu et alj|2003l) and to the "Velocity Jump" method of 
IVazza et alj d2009l) . We decided, however, to use pressure as main 
variable of the method, since the CRs pressure is also the physical 
quantity used in the code (see Sec l2.3.3l for a discussion). 

Our code first selects candidate shocked cells by requiring 
that the flow in the cell is converging: V ■ v < 0. Then the 3— 
D distribution of cells around the candidate cell is analysed with 
one-dimensional scans along the three axes, and it is checked that 
the gas temperature, T, and the gas pseudo-entropy, S — P ■ p -7 
(where p is the gas density and 7 is th e adiabatic index ) change in 
the same direction, VS • VT > (e.g. lRvu et aT]|2003l) . 

The gradient of temperature then sets the candidate post-shock 
and pre-shock cells. The Mach number is then evaluated from 
the information of the pressure jump between cells given by the 
Rankine-Hugoniot condition: 



Pi 



2-yM 2 - 7 + 1 
7+1 



(1) 



where Pi (Pi) is the gas pressure where the temperature 
across the candidate cells is the highest (the smallest). By construc- 
tion, we consider pre-shock and post-shock cells with a stencil of 3 
cells around the minimum of the 3-D divergence. 

In the case of candidate shocked cells with several Mach num- 
bers along the different axes, and in case of multiple candidate 
shocked cells along the same direction, we "clean" the preliminary 
map of Mach numbers by retaining only the maximum Mach num- 
ber provided by the algorithm, for the considered patch of cells 
(e.g. lRvu et al ] |2003h . Obliquity in shocks can also lead to a small 
additional smoothing of the shock trans ition across a few cells (e.g. 
lTaskeretal.ll2008l : ISkillman et ai]|2008l) . We verified that in the 
case of oblique shocks in 3-D, the Mach number reconstructed 
by our procedure is equal, within a ~ 5 — 10 percent even for 
strong shocks, to that measured for perpendicular shocks having 
the same input Mach number. Given our cleaning procedure on 
multiple shocked cells along the same scan direction, our code also 
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Figure 1. 2-dimensional slice showing log 10 (T)(top panels, in units of [K]) and log 10 (M) reconstructed in run-time by our code (bottom panels), for the 
full simulated cubic volume with the side of 80 Mpc at 512 3 at z = 0. The left panels show the maps for a slice of depth 300 kpc, while the right panels show 
the (volume-weighted) projected maps for a line of sight of 20 Mpc. 



makes sure than the reconstructed pattern of shocks, which is used 
for the injection of CR energy, is always 1 cell thick (see the first 
two maps of Fig[3}. 

Based on our tests dVazza et"al]|2009h . the use of a stencil of 
~ 3 cells ensures in general the best reconstruction of the shock 
jumps in cosmological runs. Also, this avoids the contamination 
of density/temperature fluctuations associated to substructures in 



clusters (e.g. small gas clumps) when the Mach number is com- 
puted. Usually, only very strong shocks (M 3> 10) are spread 
across more than 3 cells in cosmological simulations, so only in 
these cases we underestimate M slightly. However, these shocks 
are very rar e and not relevant for the production of CR energy in 
clusters (e.g.lRvu et alj2003l ; |Pfrommer et al.l2006l : ISkillman et all 
l2008tlVazzaetalJl2009h . Given the shape of the acceleration effi- 
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ciency we use here (see Sect l2.3.fl and Fig0, the exact value of 
shock strength is not crucial for M > 20 — 30, since the assumed 
efficiency function saturates for very strong shocks. 

Figure Q] shows the results of our run-time detection scheme 
for a cosmological box of 80 Mpc simulated with a fixed grid res- 
olution (Runl_h in Tab l3.3l Ax = 156 kpc). We show volume- 
weighted maps of gas temperature (top panels) and of Mach num- 
ber reconstructed at run-time by our algorithm (bottom panels) for 
a thin slice of 300 kpc (left panels) and for a much larger slice 
with the thickness of 20 M pc (right panels). S i milar to what is 
found in other papers (e.g . Miniati et al. 2001 ; iRvu et al.l [20031 ; 
Skillman et IDI2QP8; Vazza et al. 2009), the outer regions of clus- 
ters and filaments are surrounded by strong, M 3> 10, regular pat- 
terns of shocks, while the hotter innermost regions are crossed by 
more irregular and weaker shocks. Even when seen in projection 
along a large simulated volume (right panels), our reconstructed 
shocks trace sharp surfaces surrounding galaxy clusters and fila- 
ments, representing efficient and well-confined sites for the injec- 
tion of CRs inside large-scale structures. We remark that this is 
already at variance with cosmological SPH simulations, where the 
outer shocks are much smoothe r due to the large smoothing length 
outside large-scale struc tures JPfrommer et al.l 120071 ; iHoeft et all 
l2008l : IVazzaetal.ll2011bl) . 

In previous work, we demonstrated that an algorithm based 
on velocity jumps (VJ) between cells may present small advan- 
tages in the reconstruction of shock waves outside of clusters, ow- 
ing to a smaller scatter associated with velocity-bas ed measure- 
ment s compared to temperature-based measurements ( Vaz za et"al] 
2009). However, the vectorial analysis of the 3-D velocity field 
and the "cleaning" algorithms associated with the procedure, re- 
quires a substantially greater amount of calculations compared to 
scalar-based methods. This is not an issue for any post-processing 
analysis of shock waves, but it makes a significant difference in per- 
formance when the computation is performed in run-time (due to 
3-D loops and hig her memory usage). For this reason, we consid- 
ered a pressure-based method a good compromise for the run-time 
use needed here. 

In Fig|2]we show the differential volume distribution of Mach 
numbers for the full simulated cosmological volume, measured 
with the two methods. The distributions of shocks are divided ac- 
cording to the cosmic environment, depending on the total matter 
over density of shocked cells, S p = ptot/pcr (where p is the to- 
tal gas+DM density and p cr ~ 9.31 • 10~ 30 g/cm 3 is the critical 
density of the Universe: voids, 8 P < 3, filaments, 3 ^ S p < 30, 
and g alaxy clusters and their outskirts, S p ^ 30, as in I Vazza et al.l 
l2009h . 

The distributions are very similar to what has been found in 
our previous work base d on fixed resolution simulations in ENZO 
(e.g. IVazzaetalJl2011bh . At the scale of galaxy clusters outskirts 
and galaxy clusters virial volumes (S p > 30) the two methods yield 
nearly identical results. 

In runs where adaptive mesh refinement (AMR) is turned on, 
the Mach number is computed at all refinement levels, using a sten- 
cil of cells at the same level. The Mach number estimated at the 
highest resolution is always preferred in the case of several Mach 
numbers for the same cells, obtained with different AMR levels. In 
this paper, we employ adapt ive mesh refinement using the criteria 
outlined in IVazzaetalJ (2009). This means that an AMR criterion 
based on velocity jumps in 1-D is added to the usual criteria of 
gas and DM over density (e.g. O'Shea et al. 2004), ensuring that 
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Figure 2. Volume distribution of Mach numbers in a simulated cosmolog- 
ical box of side 80 Mpc (RunlJi, see Table l3.3V The solid lines report the 
results for the run-time shock finder employed in this work, the dot - dashe d 
lines are for the velocity jumps method introduced in Vazz a et alj J20O9I) . 
run over the same data. The different colours show the results for different 
cosmic environments. 

virtually all shocks with (M > 2) within the AMR regions are 
re-sampled down to the maximum available resolution. Hence, al- 
most the entire population of shocks contributing to CR accelera- 
tion (Sect l2.3.Tt is reconstructed and measured using the highest 
available resolution in the domain. Figure [3] shows the pattern of 
Mach numbers reconstructed in run-time by our method (left panel) 
for a galaxy cluster. It also shows the energy flux of CRs and the 
gas energy for a runs employing 5 levels of mesh refinement (the 
slice has a side length of ~ 5 Mpc and a depth along the line of 
sight of 32 kpc). 

In the case of a mixture of gas and CRs, the total Mach num- 
ber of the shocks can be obtained by Eq[TJ provided that the to- 
tal pressure and an effective adiabatic index (as 7 e g = (7-Pg + 
7crPcr)/(-P g + PrA lEnBlin e"tai]|200l where 7cr = 4/3). 

2.3 Method for Cosmic Rays 

In the following sections we describe our numerical implementa- 
tion of the injection, the spatial advection and the pressure feedback 
of CRs. 

In this work, we consider only CR protons, for which to 
first approximation radiative and Coulomb losses and the inelas- 
tic proton-proton collisions can be neglected, for the time-scales of 
interest here (e.g. Miniati 2007). The extension to CR electrons and 
the modelling of radiative, Coulomb losses and proton-proton colli- 
sions for both species will be investigated in future work. Also, we 
only consider CRs generated at large-scale structures shocks. The 
inclusion of other sources of CRs, such as supernovae or AGN, as 
well as the co upling of our metho d for CRs to the new MHD ver- 
sion of ENZO dCollins et alj2010f) will also be the subject of future 
work. 

In summary, in our implementation we inject the energy of 
CR particles at shocks during run-time assuming Diffusive Shock 
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Figure 3. 2-dimensional slices showing the Mach number (left, color coding as in the bottom left panel of Fig[T}, the energy injected into CR during the last 
time step (central panel, in logio [codeunits] ) and gas energy (same color coding of central panel) for a slice with side of ~ 5 Mpc along a line of sight of 32 
kpc. The simulation employs nested initial conditions and 5 levels of mesh refinement. The fact that only cells with M > 2 are used for the injection of CRs, 
explain the small difference in morphologies between the left and the central panel in the Figure (see Sec l2.3.ll for details). 
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Mach 

Figure 4. Tre nd with the Mach nu mber of the acceleration efficiency at 
shocks for the iKang & Jonesl2007t) efficiency (rj(M), in red), for the stan- 
dard thermalization efficiency in the post-shock (8(M), in blue) and for 
the ratio of the two (black). In this plot, the kinetic energy flux across the 
surface of the shock is normalized to one. 

Acceleration l |2.3.1| l, we remove this energy from the thermal gas 
energy pool l !2.3.2t . we advect it in the simulated volume along 
with the thermal gas energy l !2.3.3t , and we compute its pressure 
feedback by adding to the total pressure l |2.3.4t . In order to pre- 
serve numerical stability, new time stepping criteria are added to 
the standard ones (Sec l2.3.5t , 

2.3.1 The injection of Cosmic Rays 

In our treatment of DSA, we rely on the results of numeri- 
cal, mostly one-dimensional studies of CR acceleration at shock 
waves as a function of background plasma parameters (e.g. 



Malkov & O'C Drunfc00ll ; lKang & Jonesl2007l:lKang et alj2009t 
Caprioli et alfcOld) 

The injection of CR energy in 1-D quasi parallel shocks has 
been studied in detail with advection-diffusion equations for a va- 
riety of shocks, for velocities in the range 150 v s ^ 4500 km/s 
and for pre-shock temperatures of 10 4 — 10 6 K. In these models 
the partial dissipation of shock kinetic energy into the excitation of 
Alfven modes is accounted for, while in our treatment we neglect 
the presence of pre-existing CR energy in computing the injection 
of new CR energy. This is clearly a limitation of our treatment. 
However, self-consistent recipes to compute the estimated acceler- 
ation efficiency as a function of both the shock Mach number and 
the ratio of CR energy to gas thermal energy at every shock are 
presently unavailable, since they require the f ull self-consistent so- 
lution of diffusion-c onvection equations (e.g. IKang & Jonell2007l: 
ICaprioli et alfeoiol) . 

Based on these theoretical papers, we assume that a small frac- 
tion of the incoming protons is instantly accelerated at the dissi- 
pative gas sub-shock to a speed larger than the post-shock sound 
speed, and that it contributes to the CR pressure inside the cell. 

For every detected shock we compute the co-moving shock 
speed v s = Mc s and the resulting energy flux of accelerated CR 
protons, as a function of the Mach number: 

0cr = rj[M) ■ ~y~ = n{M) , (2) 

where p u is the co-moving up-stream gas density and r](M) 
is a function o f M, w hose numerical approximation can be found 
in IKang et alj J2007h . For simplicity we restrict the acceleration 
of CRs to M ^ 2, given that the acceleration of CRs for very 
weak shocks is exp ected to be extrem ely small, and it is still very 
poorly constrained ( K ang & Rvvj201of). T he trend with Mach num- 
ber of t?(M) in the Kang & Jonesl d2007h model is shown in Fig|4] 
along with the trend of the standard thermalization efficiency from 
Rankine-Hugoniot jump conditions (<5(Af)) and the ratio between 
the two. 

To obtain the injected energy density of CRs within each cell, 
we integrate the CR energy flux over the time step: 
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where p is gas density in the cell, dx\ is the cell size at the 
refinement level 1 and dt\ is the time step of the simulation at that 
refinement level. The CR energy density, e cr , is therefore our main 
variable to simulate the evolution and dynamics of relativistic par- 
ticles in the two-fluid model, and it represents the integrated energy 
density over the energy/momentum spectra of the CR distribution, 
e cr oc JE- N(p)dp oc Jp 4 ■ f(p) ■ (l+p 2 )-^ 2 dp. 

In every cell, the contributions of the CR energy fluxes from 
all directions are added, and the CR energy at the various AMR 
level is handled as all other primary quantities in ENZO (see 
ICollins et alj2010l for a review). 



2.3.2 Reduced thermalization at shocks 

The transfer of energy from the thermal gas to CRs causes a 
reduced thermalizati on efficiency in the post-shock region (e.g. 
lAmato & BlaslT2 005). This effect can be modelled by removing the 
corresponding energy flux of accelerated CR from the post-shock 
region, thus ensuring energy conservation where Eq[3] is applied. 

In ENZO, th e "dual energy formalism" (RyueTai] 1 19931 ; 
iBrvan et al.|[l995h is used, which means that whenever the inter- 
nal gas energy, e g , cannot be calculated correctly as the difference 
between total, etot = e g + v 2 /2 (where v is the modulus of the 
velocity field), and kinetic energy (due to the fact that etot 3> e g 
in highly supersonic flows, and that there round-off errors domi- 
nate) the total energy is evolved. Since the acceleration of CR af- 
fects the thermalization efficiency of post-shock regions, we need 
to update the two energies in a consistent way. Once e cr has been 
computed (EqO, we update the values of the total energy and of 
the gas internal energy within the cell according to: e' g = e g — e cr 
and e' tot = etot — e cr . This step is performed before starting the 
Riemann solver of the PPM method at the following time step. In 
cosmological runs with out heating from the re-ionization field (e.g. 
lHaardt & M adau 1996), the temperature of the background gas can 
reach un-realistically low te mperatures, T < 1 K, with extremely 
strong shocks (M > 1000. IVazza et al .1120091) . In these cases, one 
may end up in the situation of having e g — e cr < in some cells, 
due to similar round-off errors as in the estimate of e g . Considering 
that this problem arises only in extremely rarefied regions and in the 
absence of a re-ionization background, we just fixed it by forcing 
our version of ENZO not to inject new CRs if this is going to pro- 
duce negative values of e g or etot • This fix has no consequence on 
all thermal and non-thermal statistics of our simulated large-scale 
structures, and this problem does not arise when realistic models 
for background radiation from re-ionization are adopted. 



2.3.3 Spatial advection of CR energy 

Once injected, the CRs are spatially advected assuming they are 
"frozen" into the gas velocity field. This follows from neglecting 
the effect of the spatial diffusion of CRs during the simulation. 

Run-time processing of spatial diffusion of CRs is computa- 
tionally expensive, and requires the knowledge of the B field topol- 
ogy and of MHD modes in the simulated volume. However, accord- 
ing to quasi-linear theory, the spatial diffusion is not expected to 
play a role for the physical scales we study here, L > 25 — 100 kpc, 
given that the typical diffusion time of ~ GeV particles in the ICM 
is of the order of ~ 1 — 10 Gyr for scales of L ~ 100 kpc assuming 
B ~ 0.1 - lfiG (e.g. lBlasi etai1l200l as a recent review). Since 
the bulk of CR energy is injected at z ^ 1 (Sect |4. it , the effect of 
spatial diffusion is small compared to the typical scales analysed in 



our simulations, ~ 100 — 1000 kpc. However, the distributions of 
CRs can be slightly modified by the action of diffusion, which may 
smooth our distributions on L^g ~ 50 — 100 kpc. Since the self- 
consistent inclusion of diffusion (and co nduction) processes in th e 
simulated ICM is still challenging (e.g. iRuszkowski et al.ll20lH) . 
and since a self-consistent magnetic field is not included in our ver- 
sion of ENZO, we consider the above as an unavoidable limitation 
of our present method. 

No other dynamical exchange between CRs and thermal gas is 
assumed to take place (e.g. Coulomb losses), and the coupling be- 
tween CRs and baryonic gas proceeds via the momentum equation 
(see below). 

The energy of CRs is advected in space by solving: 



dE c 
dt 



+ V • (£crV) + P cr ■ V ■ V = 0, 



(4) 



where E CI = pe cr , and where we assume a constant pseudo- 
adiabatic index for the CR energy (e.g. Kang & Jones 1990; Math- 
ews & Guo 2011). Consequently, cosmic ray pressure, P cr , and 
energy , e cr , are related via 



-Per = (7cr - l)pe cr = (7cr — l)E cr • 

Equation[4]is numerically solved as: 



(5) 



K 

cr, 1 + 1/2 

where 



EZi + dt/dx[{E 



n+l/2 



-1/2 



+ 1/2 
-1/2 



1/2 

i+1/2 



) + Per 



+ 1/2 
-1/2 



n+l/ 2) 
J i+l/2 > 



7c: 



{Eci,i-l/2 + Ecr,i- 



1/2) 



(6) 



is the "centred" CR pressure between the the two faces of the cell 
(e.g. Bryan et al. 1995). 

Equation[6]is computed similarly to the total energy and inter- 
nal energy equations in ENZO, via one-dimensional sweeps along 
each direction and using the information of the fluxes, wJL]/ 2 and 
"1+1/2' reconstructed by the PPM method (e.g. Bryan et al. 1995). 

The effects of cosmological expansion on CRs is treated sep- 
arately, by updating the CR energy density, e cr : 



de CI 
~~dt 



(7cr - i)a 

O 6c 



(7) 



where a and d are the cosmological scale factor, and the scale factor 
rate of change at each time step (e.g. Bryan et al. 1995). 
The real value of the CR adiabatic index 



7cr = I 



d log Per 

dlogp 



(8) 



cannot be derived self-consistently, since this would require 
the knowledge of the distribution of the CR momenta, f(p). De- 
pending on the shape of the momentum distribution and of the min- 
imum momenta, g m i n , this may range from 7* r = 4/3 to 7* r = 5/3 
(e.g. Ensslin et al. 2007). In numerical simulations of DSA, where 
the diffusion-convection equation for the evolution of the f(p) is 
solved in run-time, and coupled to the dynamical equation of the 
gas evolution (e.g. Kang & Jones 2002), the adiabatic index can 
be computed self-consistently. In our case, the value of 7* r = 7 cr 
can only be assumed a-priori. In the following, we will produce 
tests and simulations where we fix the adiabatic index to its most 
extreme value of 7 cr = 4/3. 
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Figure 5. Shock tube test for P cr ,i = 6, P gi i = 10 and p\ = 9, and assuming an acceleration efficiency of r\ 0.05 at the shock front. The left panel 
shows the profile of gas density, the central one the profiles of P g , P cr and Ptot, while t he right one sh ows the profile of velocity. The diamonds represent the 
solution for our code, while the solid red lines shows the numerical solution obtained by Miniati 1 2007). 
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Figure 6. Shock tube test for the various physical modules available to our method. From left to right, we show: gas density, CR pressure, total (thick lines) 
and gas (thin lines) pressure, velocity field. The initial conditions are as in Fig|5] and r\ = 0.02 is assumed at the shock front. 
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Figure 7. Zoomed version of Fig[6]for the CR pressure at the shock front. 
The meaning of colours is as in Fig [6] 



2.3.4 Dynamical feedback of CR 

The dynamical coupling between CRs and the thermal gas is done 
by updating the gas pressure with the total pressure from the mix- 



ture of gas and CRs within every cell. At each time step, the total 
thermal gas pressure, P s = p(7 — l)e g (with 7 = 5/3) is replaced 
with the effective pressure in the cell, P e g = P g + P CI , where P CI 
is given in Eq|5] 

The effective pressure P e g is then fed to the Riemann solver in 
the next time step, and the new "left" and "right" states around each 
cells are reconstructed based on the new effective pressure. The 
momentum equation from these states is then solved in the usual 
way in ENZO, and the resulting velocity field is used to update the 
evolution of all gas quantities in the simulation, as well as of CRs 
energy field, E CI (Eq[4]l. 

In order to achieve numerical solutions without spurious oscil- 
lations, it is crucial that all fluxes associated with the total pressure, 
P ff, and with the CR pressure, P cr , fed into the Riemann solver are 
passed through the same slope limiters in order to keep the sharp- 
ness of solutions around shocks and in strong rarefaction regions 
(e.g. Colella & Woodward 1984; Bryan et al. 1995). Our tests sug- 
gests that this is particularly important when the energy density of 
CRs is of the order of e cr ~ 0.1e g , or larger. 

Furthermore, when the "left" and "right" edge states of the 
Riemann problem are computed for the solution of the fluxes, the 
appropriate wave velocity is used; this implies replacing c s = 
y/jPg/p with 



7-Pg + 7cr-Pci 



(9) 
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in the associated routines of the PPM algorithm, as suggested 
in Miniati (2007). 

2.3.5 Time stepping 

To ensure numerical stability in the presence of feedback from CRs, 
we impleme nt 2 additional const raints on the time stepping criteria 
of ENZO (see lCollinsetal.ll2010l for a review). 

First, we modify the Courant condition on the hydro- 
dynamical time step by replacing the standard speed of sound with 
c' B (as defined in Eq[9j: 

A t = m in(-^), (10) 

where v is the 1-dimensional gas velocity along one axis, and the 
minimum is taken over the stencil of cells at each refinement level. 
We need also to ensure that the internal energy correction (see 
Sec l2.3.3t . due to the reduced thermalization process, is not larger 
than a fraction of the gas energy within the cell, since this would 
cause numerical instabilities. We achieve this by imposing: 

At cr < e cr • e g • (11) 

9>cr 

which follows from the requirement that e cr e s ^ e cr at each cell. 
Our tests suggest that e cr = 0.02 — 0.04 is a suitable choice to en- 
sure that the modification of the thermal energy in the post-shock is 
not dramatic during each time step, also in the case of the complex 
3-D flows (which may cause fluxes of accelerated CRs across sev- 
eral faces of a cell). Compared to the standard cosmological sim- 
ulations at fixed grid resolution, the use of these additional time 
stepping criteria usually produce an increase of ~ 15 percent in the 
total CPU time by the end of runs. 



3 TESTS IN ONE DIMENSION. 

3.1 Shock Tube with pre-existing Cosmic Rays 

We performed one-dimensional shock tube tests assuming a mix- 
ture of CRs and thermal gas in the initial setup, and investigating 
the behaviour of our scheme at different Mach numbers. 

The parameters of the first shock tube test are identical to 
those reported in Miniati (2007): the initial left state is filled with 
a mixture of gas and CRs, with pressure (in code units) P gi i = 10, 
Pci,\ = 6, and density pi — 9, while all the fields on the right state 
are set to 1 (7 cr = 4/3 is assumed everywhere for CRs). 

In Fig[5]we show the comparison of our test with the numeri- 
cal solution of Miniati (2007), at approximately the same epoch. In 
this case, all physical modules in our version of ENZO with Cosmic 
Rays are activated: advection, injection, pressure feedback and re- 
duced thermalization at the post-shock. The agreement of our run 
with the results of Miniati (2007) is excellent, with a constant total 
pressure level in the post-shock region. The contact discontinuity 
at the rarefaction wave is reconstructed in a few cells (which is 
common in the PPM method), while the transition of the shock is 
reconstructed in ~ 2 cells. 

The numerical implementation of our CR code also allows us 
to investigate the role of each physical module separately. We show 
the effect of each module in FigHj] where the same initial condi- 
tions of the previous shock tube were re-simulated by modelling: 
only the advection of CRs is modelled (run SO); the advection of 
CRs and their pressure feedback onto the gas component (SI); the 
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Figure 8. Behaviour of the gas pressure (dots), CR pressure (diamonds) and 
total pressure (solid lines) for shock tube tests leading to a M 1.5 (blue), 
M fS 3 (green) and M fs 5 (red) shocks. Only the values of CR pressure 
are shown for clarity; the x-axis is in units of the total box length. 

injection of new CR energy at the shock (52); the whole set of CR 
modules (S3). 

When we allow CRs to have pressure feedback on the gas (SI), 
the total pressure of the system is increased and the behaviour of 
gas density, velocity and of the total pressure are scaled-up ver- 
sions of those in 50. However, the behaviour of gas pressure at the 
shock front is modified by the presence of the CR pressure. The 
injection of CR energy at the shock front in 52 (here r\ ~ 0.03) 
slightly increases the total pressure at the shock compared to SI. 
When including the reduced thermalization (S3), the total pressure 
at the shock front is smaller compared to 52, and as a result the 
speed of the shock front is decreased. Figure|7]shows a close-up of 
the behaviour of CR pressure at the shock front. Only very small 
oscillations (at the percent level) are present in the post-shock so- 
lution of our CR fluid. 

In all runs where the CRs are allowed to exert a pressure on 
the gas (57, 52 and S3), the solution in the rarefaction region does 
not vary. This confirms that the treatment of CR acceleration at the 
shock does not affect any other region in the shock tube. 

In a second set of shock tube tests, we varied the initial con- 
ditions by increasing the gas pressure on the left side in order to 
produce stronger shocks. In Fig[8]we show the results for shock 
tubes leading to M f» 1.5, 3 and 5 and assuming r\ ~ 0.3 to high- 
light the effect of CR injection at the shock. Even for the strongest 
shock, the post-shock total pressure shows no significant spurious 
oscillation, and the shock edge is still reconstructed in 2-3 cells at 
most. 

3.2 Shock Tube without pre-existing Cosmic Rays 

Next, we focus on the injection of CR energy at the shock front 
by simulating similar 1-D shock tubes as in the previous case, but 
assuming no pre-existing CR energy at the beginning of the simula- 
tion. Also, we increased the values of initial gas pressure on the left 
side of the box (P g ,i =10, 20, 30, 60 and 120). In these tests, we let 
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Figure 9. Behaviour of the gas pressure (dots), CR pressure (diamonds) and 
total pressure (solid lines) for several shock tube tests without pre-existing 
CR energy. Only the cell-wise values for the total pressure are shown for 
clarity. 
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Figure 11. Comparison of the total pressure of the Zeldovich collapse test 
in Z3 employing fixed resolution and of the 73 with 4 levels of AMR. The 
different lines show: gas pressure (solid), CR pressure (dot-dashed) and the 
total pressure (dotted). 



rj(M) vary with the shock Mach number according to the model of 
Kang & Jones (2007), which we will use in our cosmological runs. 

The results are shown in Fig(9] At the strongest shocks, the 
freshly injected CR pressure becomes comparable to the post-shock 
gas pressure (P CT ~ 0.3P g ). The total pressure behind the shock 
still remains constant, despite a small amount of oscillations of nu- 
merical nature. Again, the shock edge is reconstructed within 2-3 
cells at most in all tests. 

To avoid any smoothing of the shock region, which may cause 
spurious effects when the non-linearity of CR acceleration schemes 
is coupled to the system, Miniati (2007) proposed the application 
of a hybrid Godunov-Glimm's scheme to ensure the sharpest re- 
construction of shocks in the simulation. 

In order to limit the effect of a broadened shock profile on the 
estimate of the injected CR energy in all our tests here and in the 
production runs (Sec[4j, our run-time shock finder (Sec l2.2t mea- 
sures Mach numbers in the simulations across 3 cells. This was mo- 
tivated by our previous work on cosmological shock waves in ENZO 
simulations (Vazza, Brunetti & Gheller 2009), where we studied the 
best stencil of cells needed to correctly capture shock waves in the 
PPM version of ENZO. With this approach, the effect of numerical 
smearing of reconstructed shocks only has a very limited impact on 
the estimated injection efficiency of CRs, since M is always mea- 
sured with a stencil of cells close to the (small) numerical smearing 
of shocks in the PPM method. 

3.3 Zeldovich Pancake 

In order to test our implementation of cosmic ray physics in a sim- 
ple cosmological framework, we run a set of Zeldovich Pancake 
test and compare the results with and without the dynamical role 
of CR energy injected at shocks. The initial conditions (z\ = 50 
in our case) for this test assume a purely baryonic Universe, with 
a uniform initial pressure and density, and a converging velocity 
field with v — at the centre of the domain. With this setup we 



simulate the evolution of a co-moving box of Lbox = 20 Mpc/h 
(h — 0.72) with a grid of 256 cells. No CR energy is present in 
the box at the start of the simulation. Here results are presented for 
the simple collapse test (run Zl), and for the full set of physical 
modul es for CRs (Z3) . In all cases the acceleration efficiency is the 
one o f iKang & Jonesl (2007). Additionally, in order to test the reli- 
ability of our method with adaptive mesh refinement, an AMR run 
with a 16 times higher spatial resolution has been performed (run 
Z3+AMR). 

The panels in Figure [10] show the evolution of gas density, gas 
pressure, CR pressure, total pressure and total energy in runs Zl 
and Z3. By the end of the runs the profile of the "pancake" is very 
similar in all cases, with sizable differences at the percent level 
which are evident for z < 2, after enough CR energy has been 
injected by the expanding shocks. The gas density and gas pressure 
of Z3 are slightly smaller inside the structure and slightly larger 
in the external regions (at ~ —0.8 — 1.2 Mpc/h), while the total 
pressure profile is very similar in the two cases, with the exception 
of the outer region where run Z3 shows an excess. The reason for 
this behaviour is that, when output at equal times is compared in 
runs Zl and Z3, the outer shocks in the simulation with CRs are 
found to have expanded to a slightly larger distance, due to their 
enhanced total pressure. The combination of the softer equation of 
state of the gas+CR composite adopted in Z3 at the outer regions, 
and the total larger pressure jumps felt by the gas accreted onto 
the pancake makes this mechanism very stable, since both effects 
concur to have a more efficient injection of CRs in the system, due 
to the increase of the Mach number of external shocks. We will 
discuss this mechanism again in Sec.4 and in the Appendix, with 
the results of fully cosmological runs in 3 dimensions. 

By the end of this run, the pressure ratio between CRs and gas 
is Pcr/Pg ~ 0.6 at the external regions and P CT /P S ~ 0.05 — 0.1 
in the innermost region of the pancake. 

In FigQjjwe show the effect of AMR (4 levels of refinement, 
triggered by the local gas overdensity) on the same collapse test. 
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Figure 10. Time evolution of the profiles of gas density, gas and CR pressure (dotted lines in the second panel), effective adiabatic index (bottom left) and 
total pressure (bottom right) for the Zeldovich collapse test run with no CR physics (solid lines, run Zl) and with CR physics (dot-dash, ran Z3). 



The solution recovered by the end of the run is quite similar to Z3, 
provided that the larger spatial resolution here enables a sharper 
modelling of the outer shock, and leads to a larger injection of 
CRs. These two effects modify the shock speed in a non-linear way, 
making the synchronization of the two runs not trivial. The energy 
ratio measured at the outer region is larger compared to Z3, yield- 
ing P CI /P S ~ 0.5. Based on the profile of P g and P cr along the 
pancake, we conclude that the previous results at coarser resolution 
are stable against such a dramatic increase of resolution, even if 
the detailed distribution of P cr at the shock edge may slightly vary 
with resolution. As in the previous case, this difference can be par- 
tially ascribed to small timing issues in the comparison of the two 
runs, and partially to the slightly enhanced amount of CRs injection 
measured at higher resolution. 



4 RESULTS IN COSMOLOGICAL SIMULATIONS 

We produced a large set of cosmological simulations in order to 
study the evolution of CRs accelerated at accretion shocks in large- 
scale structures, and their dynamical feedback on the clustering 
properties of the thermal baryon gas. In a first set of simulations, 



we studied the large-scale properties of CRs using fixed grid sim- 
ulations (Sect |4. it , while in a second set we studied the details of 
the distribution of CR energy in galaxy clusters re-simulated at high 
spatial resolution with AMR (Sect |4.2| l. 



4.1 Runs with fixed resolution 

In a first set of cosmological runs we simulate the evolution of a 
box with a co-moving side of 80 Mpc, in a "concordance" cos- 
mological ACDM model, with fi dm = 0.226, n h = 0.044 and 
SIa = 0.73. The normalization of the primordial index of density 
fluctuations was set to erg = 0.8 and h = 0.72. The initial red- 
shift of the simulations is z- m = 30. All runs are non-radiative, and 
neglect all non-thermal phenomena other than CR physics, such 
as energy release from stars and AGN or magnetic fields. When 
re-ionization in the early Universe is modell ed, we follow the run- 
time recipe introduced in lVazza et al.l (20101. where a simple pre- 
scription for a temperature floor in the simulation is adopted in the 
redshift range 7 ^= z ^ 2 in order to mimic the re-ionization, as in 
lHaardt & Madaul dl996t) . Despite the moderate resolution achieved 
in these runs (Ax ^ 122 kpc/h), they allow us to test the role of 
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Figure 12. Comparison for the distributions of Mach number (upper left panel), gas temperature (upper right), gas density (lower left) and P cr /-Pg as a 
function of gas over-density (lower right) for runs Runl, Run2, Run3, Run3_l, Run3_m and Run3_h. 
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Figure 13. Maps of CR pressure as a function of resolution, for a slice of 
20 X 60 Mpc and width 448 kpc/h, centred on the same massive galaxy 
cluster in the simulated volume (Run3, Run3_m and Run3_h are shown). 



the parameters of our models. The summary of the 15 runs at fixed 
grid resolution is listed in Tab. 13.31 

Our fiducial model for simulating the evolution of CR energy 
in our ENZO cosmological runs is Run3 (and Run3_l, Run3_m and 
Run3_h at other resolutions, see below). In this model we adopt the 
full set of CR-modules of our implementation (injection at shocks, 
spatial advection, pressure feedback and the reduced thermaliza- 
tion). The injection of CRs at shocks starts from z CI = 6, and 
for p/pcr.b 0.1 (where p cr ,b is the critical baryon density of 
the Universe). The a cceleration efficiency depends on M, as in 
iKang & Jonel d2007t) . 

In other runs we varied various parameters in order to study: 

• the dynamical role of CRs in the evolution of large-scale struc- 
ture with a simple non-radiative run (Runl). In this run, we inject 
and passively advect CRs in the simulated volume, but CRs do not 
couple to the baryoicn gas (therefore they do not exert pressure nor 
do they cause reduced thermalization in the post-shock); 

• the role of the reduced thermalization efficiency, by compar- 
ing Run3 with a run with CRs injection and pressure feedback but 
no reduced thermalization at shocks: Run2 (in this model the total 
energy at the shock is not conserved, since "new" pressure from 
accelerated CRs is just added to the system); 

• the role of spatial resolution and DM mass resolution by re- 
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Table 1. Numerical and cosmological parameters for the runs employing 
fixed mesh resolution. The columns lists: ID of the run; mass resolution for 
DM; gas spatial resolution; number of grid cells. The last column reports 
mnemonics explaining the physical setups used in the run. "A"=advection 
of CR; "/ "=injection of CR at shocks; "7"=reduced thermalization at the 
post-shock; "P"=pressure feedback of CRs. "z CT = 1, 29" mean that the 
injection in the runs was started at z = 1 and z = 29, respectively (or 
at z = 6 elsewhere), "no-re." means that the run does not model a re- 
ionization background (otherwise considered), "rj = 0.01" and "r) = 0.01" 
mean that the acceleration efficiency in the runs was kept fixed at these two 
values, independent on the measured Mach number (the r)(M) efficiency 
of Kang & Jones (2007) was used elsewhere). "M cr " means that the Mach 
number used to compute r)(M) is the one from the total gas+CRs pressure 
jump, instead of the gas pressure jump only. "p cr = 10 -3 , 10" means that 
we assumed minimum overdensity of 10~ 3 or 10 the background critical 
gas density to inject CRs (p cr =0.1 elsewhere). 
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Run2 
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5.1 


■ 10 9 


448 


128 3 


AIPT, z cr = 1 


Run3.z29 


.5.1 


■ 10 9 


448 


128 3 


AIPT, z CT = 29 
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AIPT, p cr = 10~ 3 


Run3.dl0 


5.1 
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128 3 


AIPT, p cr = 10 


Run3_mcr 


.5.1 


■ 10 9 


448 


128 3 


AIPT, M cr , Tj = 0.1 



Table 2. Main characteristics of the simulated clusters. Column 1 : identifi- 
cation number. Columns 2: total mass (A7tot = A-7dm + A7 gas ) inside 
the virial radius. Columns 3: virial radius, R v 0.77?200- Column 4: dy- 
namical state at z = 0: MM=major merger cluster; ME=merging cluster, 
RE=relaxed cluster (see Sect l4.2l for further details). 

ID A7 TO t[10 14 A7 q ] iJ v [kpc] dynamical state 



El 11.20 2670 MM 

E7 6.52 2194 ME 

E25 6.55 2369 MM 

HI 3.10 1890 RE 

H3 2.95 1710 ME 

H5 2.41 1703 MM 

H7 2.14 1410 RE 



simulating Run3 with a coarser (Run3_l) and with two finer resolu- 
tions (Run3_m and Run3_h). For comparison, we also re-simulated 
the same setup of Runl with 512 3 (RunlJi); 

• the role of re-ionization by re-simulating Run3, for the ex- 
treme assumption of no re-ionization at all (Run3_nore); 

• the effect of assuming different epochs for the start of the in- 
jection of CRs at shocks, by investigating also the cases of z CI = 1 
(Run3.zl) and z CI = 29 (Run3_z29); 

• the effects of assuming a different minimum gas density for 
the injection of CRs at shocks, studying the case of p/p cr ,b = 10~ 3 
(Run3_d0001) and p/p cr ,b = 10 (Run3_dl0); 

• the effects of assuming fixed efficiencies for the accelera- 
tion of CR at all shocks: we investigated the "toy" models where 
the fixed efficiencies of r\ = 0.1 (Run3_eta01) and r\ = 0.01 
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Figure 14. Average pressure ratio between CRs and thermal gas as function 
of the gas over-density in our runs at fixed grid resolution for 2 = 1 (top) 
and 2 = (bottom). 



(Run3_eta001) are adopted at all shocks (with no dependence on 

M); 

• the differences caused by computing the Mach number based 
on the tot al pressure (Run3_mc r) instead of gas pressure, and ap- 
plying the lKang & Jones! J2007l) efficiency. 

4.7.7 Large-scale distributions of thermal gas and cosmic rays 

We first test the CR1 modules by comparing the large-scale dis- 
tributions of Mach number, gas density, gas temperature and CRs 
pressure in Runl, Run2 and Run3 at z=0.5 (Fig|12b. The presence 
of CRs has no strong dynamical effect in the evolution of thermal 
gas on these large scales. Some differences are found at large over- 
densities (where CR feedback reduces by a ~ 5 — 10 percent the 
largest over-densities in the box), and in the temperature distribu- 
tions inside accretion shocks, T > 10 5 K. Inside large-scale struc- 
tures, Run3 presents a small deficit of gas temperature compared to 
the standard simulation of Runl, due to the fact that part of the ther- 
mal energy at accretion shocks is channelled into energy of CRs. 
On the other hand, Run2 shows a significant excess of temperatures 
within the same regions compared to Runl. The obvious explana- 
tion for this is that the total energy at shocks is not conserved in 
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Figure 15. From top to bottom: 3-D volume rendering of gas den- 
sity (in [p/Pcrb]). Pg anQ Per (arbitrary code units) for (80Mpc) 3 
in Run3_h at 2 = 0. Rendering performed with VISIT 2.3.1 

(https://wci. Unl.gov/codes/visit). 
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Figure 16. Scaling laws for the 40 most massive halos in Runlji (black 
dots) and Run3_h (blue). Top panel: scaling between the total cluster mass, 
M200 (in Mq and the mass-weighted temperature, T200- Central panel: 
scaling between M200 anQ the total X-ray bolometric luminosity inside 
R200, ^x,200- The straight lines show the best-fit for each sample. Bot- 
tom panel: average pressure ratio between CRs and gas for all clusters in 
Run3.h. 



this run, but new energy is added to the system when new CRs are 
injected. This makes outer accretion shocks stronger over time in 
Run2, and therefore more efficient in thermalizing the medium. 

An important proxy to understand the dynamical role of ac- 
celerated CRs is the pressure ratio between CRs and thermal gas, 
Pcr/Pg, as a function of cosmic over-density (last panel of Fig] 121 
The trend of this ratio is broadly similar in all runs, with a max- 
imum at the over-densities close to the critical one. In our fidu- 
cial model (Run3) this maximum is P cr /P g ~ 0.2, and smoothly 
declines to P C r/P g ~ 0.04 - 0.08 for p/p cr ,b ~ 100. This ra- 
tio is larger if the reduced thermalization at shocks is not included 
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(Pcr/Pg ~ 0.5 at the maximum in Run2) and even larger if the ac- 
celerated CRs are just passively adverted in the simulation (Runl). 

We test the effect of spatial resolution and DM mass resolution 
by re-simulating the same setup of Run3 with at a smaller (Run3_l, 
using 64 3 cells and mdm = 4.0 ■ 10 10 'Mq /h) and at two larger 
resolutions (Run3_m with 256 3 and md m = 5.1 ■ 10 9 Mq /h, and 
Run3_h with 512 3 and m dm = 8 • W 7 M Q /h). The final maps 
of CR pressure for the most massive galaxy cluster in the box is 
shown in Fig[T3] and suggests that the CR pressure at the largest 
scales in the simulation are very similarly reconstructed at all res- 
olutions. The increase in resolution, however, produces a more 
sub-structured distribution of CR pressure inside clusters. As ex- 
pected from previous studies tevu et ai1l2003l ; Pvazza et al.ll2009l 
the increase in the spatial resolution leads to a sharper re- 
construction of shock waves, and to a lowering of the Mach num- 
ber of strong accretion shocks. When averaged over the whole cos- 
mological volume, very good convergence in the volume distribu- 
tion of simulated large-scale shocks is achieved at the resolution 
of Ax ~ 200 kpc/h, corresponding to the resolution achieved in 
Run3 (Figl 12b- The average pressure ratio P CT /P g shows a very 
regular behaviour across resolutions, and a nearly converged trend 
with over-density for p/p cr ,b > 1 — 10. At the maximum resolution 
here (Aa; = 112 kpc/h) the average pressure ratio at the scale of 
the innermost cluster regions is P cr /P g ~ 0.05. In Sec l4.2l we will 
show that this trend is confirmed even at larger resolution, using 
adaptive mesh refinement for a subset of massive galaxy clusters. 

We finally tested the various assumptions related to the in- 
jection of CRs by modifying the setup of Run3. For instance, a 
suitable minimum magnetic field is required for DSA to work, and 
therefore it is reasonable that the acceleration of CRs takes place 
only at over-densities large enough to have had the opportunity to 
generate sufficiently high magnetic fields. Assuming a minimum 
over-density for the injection of CRs in the simulation also reduces 
the number of un-necessary computations. 

Figure [T4l summarizes the differences in the final distribution 
of Pcr/Pg. as a function of gas over-density in the simulated vol- 
ume, for all runs with varied parameters of CR acceleration but 
identical spatial and DM mass resolution, at z — 1 and 2 = 0. 

The differences are somewhat larger at z = 1 (top panel) and 
become less significant for large-scale structures at z = (bottom). 
At all epochs, the maximum of the pressure ratios is always found 
at p ~ Pcr.b, meaning that the accretion regions of large-scale 
structures are the locations where CRs can have the largest dynam- 
ical role. As seen above, the pressure ratio inside large-scale struc- 
tures is usually much smaller, P cr /P g ~ 0.1 for p/p cr ,b ^ 10 2 , in 
all tested models. 

When the acceleration efficiency depends on the Mach num- 
ber, we find that the pressure ratio for all over-densities p/p C r,b 1 
is almost independent of the assumptions (on the initial epoch, min- 
imum over-density for injection, re-ionization, etc.). Differences of 
one order of magnitude in the pressure ratio between the differ- 
ent re-simulations can be found only for densities below the criti- 
cal one. There the effects of re-ionization and of the assumed first 
epoch of DSA are still present at z — 0. 

In the runs using a fixed acceleration efficiency the two trends 
for p/p C T,b > 1 are rescaled version of the same profile. In partic- 
ular the model with the fixed efficiency of r; = 0.1 lies j ust be low 
our fiducial model with the efficiency o f iKang & J ones (2007j) for 



r/(M). This suggests that at z — the bulk of the injection of en- 
ergy in CRs happens with an average Mach number of M ~ 2 — 3 
(typical of internal merger shocks), which are characterized by an 
acceleration efficiency of r\ ~ 0.1 — 0.2 (see also Fig@. This is 
also in line with our previous estimates based on post-processing 
the time evol ution of the energy processes at cosmological shocks 
in ENZO runs dVazza et alj2009h . 

Given the present large theoretical uncertainties conditions of 
DSA in the early universe, and in presence of small gas density and 
very weak magnetic fields, it is reassuring that such uncertainties 
contribute only modestly to the final distribution of CR energy in- 
side all structures, for z < 1. In essence, we can assume that the 
average distribution of CR energy inside large-scale structures is 
almost totally dominated by the dynamics of accretion shocks and 
merger shocks in the hot and dense baryon plasma, within the last 
- 8 Gyr. 

This ensures that, at least to first approximation and limited to 
regions of radii ~ 2 — 3_R v ir around galaxy clusters, the statistics 
provided by our method for CRs is very robust, once an accelera- 
tion scenario for CRs is specified. 

4.1.2 Scaling relations for clusters 

We first studied the global effects of CR dynamical feedback on a 
statistical sample of galaxy clusters, extracted from the two most 
resolved runs at fixed resolution (Runl_h and Run3_h). 

The spatial resolution of 1 12 kpc/h (which also represents the 
softening for the computation of the gravitational force) is suitable 
to study the outer accretion regions of galaxy clusters and groups, 
but not to the resolve in detail the mixing in the innermost cluster 
regions, and the evolution of cluster cores. For this reason, we con- 
sider the sample of clusters useful to study the effect of CR dynam- 
ics in the region inside -R200 (where -R200 is the radius enclosing an 
average density of 200 times the cosmic critical density). Clusters 
were extracted with a hal o finder based on spherical over-density 
(as in iGheller et aljfl998h . A more detailed study of the distribu- 
tions of CRs inside clusters, at a higher resolution, is presented in 
Secl42l 

Our volume of (80Mpc) 3 contains ~ 50 halos with a total 
mass larger than 1O 13 M at 2 = (8 with M > 1O 14 M ). The 
panels in Figure [15] give the visual impression of the number of 
halos in the simulated volume, and the large-scale distributions of 
gas and CR pressure associated with them. Since the injection of 
CRs occurs only for a large enough overdensity (p ^ 0.1p cr ,b), the 
3-D distribution of CR pressure has a much lower volume filling 
factor than that of gas pressure. 

The comparison of the T200 vs A/200 relation and of I/x,200 
vs M200 (i.e. average temperature, total mass and bolometric X-ray 
luminosity inside R200) for the most massive clusters of both runs 
is shown in the top panels of FigfTol 

Given the limitations of physical processes included in our 
runs (no radiative cooling, star formation and feedback from su- 
pernovae and AGN), the comparison of these two runs is helpful to 
highlight the effects of simulated CRs in clusters, while they do not 
have a strong prediction power on observed scaling relations. 

In the two runs the self-similarity of cluster relations is pre- 
served, but the dynamical feedback of CR energy injected at shocks 
modifies the normalization of both scaling by ~ 20 percent for 
(M2qo,T2oo) and for (M200, £-x,20o)- As will be shown in the fol- 
lowing (Sec |4.2t this is a result of a slightly modified distribution 
of gas matter density inside halos. The last panel of Figl 161 shows 
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Figure 17. Top panels: maps of gas density for slices of depth 25 kpc/h through the centre of re-simulated clusters of various masses at z - 
p/Pci b)- Bottom panels: maps of CRs to gas pressure ratio for the same regions. Each panel has a side of ss 10 Mpc. 
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the pressure ratio between CRs and gas within -R200 for all ob- 
jects of Run3_h. The majority of clusters has a ratio of the order of 
Pcr/Pg ~ 0.1, with an almost constant average ratio across one or- 
der of magnitude in mass (even though the distribution has a quite 
large scatter, and no simple trend with total mass can be detected). 



4.2 Adaptive mesh refinement runs 

A second set of cosmological simulations was produced with adap- 
tive mesh refinement, achieving almost uniform high spatial reso- 
lution within cubic boxes of side length ~ 4 — 6-R200 for several 
galaxy clusters in the mass range of 2 • 10 14 ^ M/Mq/H ^ 
2 • 10 15 . These runs are computationally fairly expensive (e.g. 
~ 2 ■ 10 4 hours of CPU time for the most massive clusters), and 
in this work we limit ourselves to the analysis of 7 re-simulated 
objects, selected from catalogues of clusters already presented in 
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Figure 18. Maps of gas pressure (first row, shown is log 1 g(P g )), CR pressure (central row, shown is log 10 (P cr )) and ratio of the two (last row. 
log 10 (P cr /P g )) for the major merger cluster H7 at four epochs. Each slice has depth 25 kpc/h and side 10 Mpc/h. 
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Figure 20. Cumulative gas mass distribution as a function of P cr /Pg, for 
clusters HI, H3 and H5 at z=0. The solid lines are for the volume inside 
^5 2P200 > the dot-dashed lines are for ^ P200 ■ 



IVazza et al| bOld) and I Vazza et all feOldh . in order to sample dif- 
ferent masses and dynamical state of the ICM. 

The use of AMR is necessary here to model the turbulent mo- 
tions in mergers, which may affect the properties of transport of CR 
energy within clusters over time. In addition, the increased resolu- 
tion in capturing shocks can enhance our modelling of CR injection 
at internal shocks. 

All objects have been re-simulated at high spatial and DM 
mass resolution starting from parent simul ations at lower r esolu- 
tion, using nested initial conditions (e.g. lAbel et all 1998). The 
mass resolution for DM inside the high resolution region where 
the clusters form is md m = 6.7 • 10 7 Mq, the coarsest spatial res- 
olution inside this region is 200 kpc/h while the peak one is 25 
kpc/h. In general a fraction of ~ 60 — 80 percent of the total vol- 
ume inside ^ 27?200 is simulated at the highest available resolution 



by the end of the simulation JVazza et alj|201fj|) . In these runs we 
triggere d mesh refinement b y standard gas/DM over-d ensity crite- 
na (e.g. lO'Shea et alj|2004l) and velocity jumps (as in I Vazza et al.l 
l2009h . 

Table [2] summarizes the main properties of each object. The 
last column indica tes the dynamical state of the cluster. Following 
I Vazza etaljj2010h we estimated the dynamical state of each cluster 
at z=0 from its total matter accretion history for z < 1.0 (in the case 
of the lowest masses systems, < 3 • 1O 14 M0, we consider instead 
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Figure 19. Left panel:radial profiles of average CRs to gas pressure ratio for the sample of re-simulated clusters with AMR at z = (solid lines), and pressure 
ratio only for the CRs injected during the last At ~ 0.1 Gyr (long dashed lines). Right panel: average effective (pressure weighted) adiabatic index for the 
same clusters at 2 = 0. 



2 < 0.6 for the same analysis, given their shorter dynamical times). 
If a system experienced an increase of total matter by 30 percent 
of more within 1 Gyr of time, we classify it as a "post merger" 
systems, or "relaxed" otherwise. In the case of systems without a 
major merger in the past, but having an ongoing merger process 
(estimated from the ratio between total kinetic and thermal energy 
within -R200) the system is classified as "merging". 

Figure [771 shows cuts through the centre of mass of each ob- 
ject of the sample, reporting the gas density (left panels) and the 
pressure ratio within each cell of the 25 kpc/h slice (right). Sub- 
stantial morphological differences are present, related to the dy- 
namical state of the cluster: in merging systems (e.g. E7, H3) the 
presence of close companions and ongoing large-scale accretions 
drive streams of enriched CRs gas in the ICM, even close to the 
main cluster centre. In post-merger systems (e.g. El, H5) the re- 
gions of large P cr /P g ratio are well associated with merger shock 
waves, which expand out of the innermost cluster regions. 

The long time evolution for one of these systems (H7) is 
shown in Figure[T8lfor the epochs of z = 0.91, z = 0.52, z = 0.23 
and z = 0.01. Two powerful merger shocks are launched just after 
the collision of the cluster cores at z — 0.91, and become powerful 
sites for the injection of new CRs in the ICM. While these shocks 
expand into the outer cluster regions, their Mach number increases, 
lust downstream of the shocks, the local pressure ratio becomes of 
the order of P CT /P g ~ 0.3, while inside the hot core region of the 
merger remnant the pressure ratio is small, P CI /P S ~ 0.02 — 0.05. 
This happens because internal merger shocks in their initial phase 
are weak, M ~ 2, and thus more efficient in heating the ICM 
than in injecting new CRs. On the other hand while they expand 
their Mach number increases due the drop in temperature of the 
ICM, and the efficiency of injection of CRs is larger. At later times, 
the cluster is characterized by a pretty regular density structure for 
< -R200, and regions of P CI /P g > 0.1 can be found only near 
the outer accretion regions, and also along filaments of cold gas 
(T ~ 10 J K) in some sectors of the cluster. 



4.2. 1 Radial profiles of CRs 

Despite the different dynamical states present in our cluster sample, 
the radial profiles of the pressure ratio within each system at z — 
are rather similar (left panel of Fig|19l> Q. They show a behaviour 
similar to that of thermal pressure, and are slightly flatter approach- 
ing ~ i?200- Within the core region of our systems, the pressure of 
CRs is ~ 0.02 — 0.04P g , and it increases to ~ 0.1P g approaching 
1 — 2 -R200 • In the same panel, we also show the pressure ratio of the 
CRs injected at shock waves during the last time step of each run 
(At ~ 0.1 Gyr). Only in the case of the merging system H3, a large 
amount of freshly injected CR pressure at a peripheral shock is 
comparable with the thermal pressure, while otherwise the pressure 
injected at shocks at late time is on average just a small fraction of 
the gas and CR pressure. The right panel of Fig[T9lshows the radial 
profile of the pressure-weighted average effective adiabatic index, 
7 e ff. Given the small ratio of CR to gas pressure in the innermost 
regions of clusters, 7 e ff is very close to the non-relativistic value of 
5/3 until the outer accretion regio ns are reached. These t rends are 
quite similar to what is reported in lPfrommer et alj d2007l) for non- 
radiative simulations with GADGET2. However, we notice that 
our profiles of CR pressure are significantly flatter, and our runs 
tend to concentrate less CRs in the innermost cl uster regions (e.g. 
IPfrommer et al.ll200l|Pi"nzke & Pfrommerl2010l) . In Sec[5]we will 
discuss the possible reasons for this. 

To highlight any dependence on the dynamical state of clus- 
ters we will focus in the following on three representative objects 
in the sample: HI (relaxed), H3 (merging) and H5 (post-merger). 
These three clusters have a similar final total mass but very different 
dynamical histories. 

For each of them we computed the cumulative gas mass dis- 
tribution as a function of Pj P g (Fig|20t. In all cases, roughly ~ 90 



We note that characteriz ing the exac t centre of the most perturbed sys- 
tems is not trivial, see also IVazz a et al. ( 2 011al) for a similar discussion in 
the case of radial profiles of turbulent energy. 
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percent of the gas mass inside R200 is characterized by a small pres- 
sure in CRs, P C r/Pg ^ 0.1. The post-merger system H5 is charac- 
terized by the largest gas mass with greater CRs pressure, followed 
by the merging system H3. The differences between the dynami- 
cal states become more prominent by considering the distributions 
inside < 27?2oo, highlighting the role of large-scale patterns of ac- 
cretion even for clusters with quite similar final mass. 

It is interesting to compare the runs with injection and feed- 
back from CRs to the standard non-radiative re-simulations of the 
same systems. Some very general trends emerge by computing the 
3-D radial profiles of gas density, gas temperature and gas entropy 
(Fig|21t: the innermost gas density is lower by a few percent in runs 
with CR physics. The gas entropy is correspondingly higher, while 
the gas temperature is slightly lower at all radii. Also the DM den- 
sity in the innermost radial bin is decreased by a similar amount. 
As shown also with fixed grid simulation (Sect H, 1,21 this is a very 
general outcome of our runs, and very similar results are also found 
for all the other re-simulated clusters of our sample. 

The reason for this small differences with standard runs is non- 
trivial. Starting from the early formation of cosmic structures, the 
injection of CRs is efficient in the outer accretion regions, grad- 
ually leading to a drop of the effective adiabatic index in these re- 
gions. This enhances the gas density in the outer regions, and leaves 
slightly less gas matter available for the inflow towards the central 
regions of clusters. While expanding into the rarefying cosmic en- 
vironment outside of structures, the relative pressure ratio of CRs 
is further increased because of their softer equation of state (in ad- 
dition, the outer thermal pressure is decreased by our modelling of 
the reduced fhermalization efficiency at the post-shock). This pro- 



cess leads overall to an enhanced total (gas+CR) pressure of the 
outer accretion regions of forming structures, compared to stan- 
dard runs without feedback from CRs. At the same physical time, 
the structures simulated with CR feedback have a more extended 
and dense envelope of accretion regions. This effect is similar to 
what was already found for the Zeldovich collapse test (Sect |3.3l l. 
In fully cosmological simulations this effect is amplified by the fact 
that the expansion in 3-D of the outer cluster shells makes the rela- 
tive increase of CR pressure more significant. We further explored 
the reasons of the above differences between simulations with and 
without CR feedback in the Appendix A. 



4.2.2 Properties of X-ray emission and and SZ signal 

Differences in the 3-D distribution of thermal gas between runs 
with and without CRs lead also to significant diffe rences in X-ray 
emission and SZ effect (e.g. lPfrommer et aU2007l) . We investigate 
here the projected emission maps of the simple bolometric X-ray 
luminosity (Lx oc J n^T^ 2 dV, where n c is the thermal elec- 
tron density and T a is the thermal electron temperature, which we 
assume equal to the gas temperature for simplicity) and of the inte- 
grated thermal Sunyaev-Zeldovich effect (Y oc J n T dV) along 
different projections of our clusters. 

In Fig|22]we show the 3 maps of Y for the merging cluster H3 
at z=0 (simulated without CR feedback, top, and with CR feedback, 
central panel), and the pixel by pixel ratio obtained by dividing 
the maps on the top panels with the corresponding ones from a re- 
simulations of H3 without CR feedback (bottom panels). 

The large-scale morphology of the SZ effect in the two cases 
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is extremely similar, with minor differences at small scales that can 
be partially explained by the different time-stepping of the two runs 
(which affect the exact 3-D position of accreted gas clumps in the 
main cluster atmosphere). The most prominent difference is related 
to the outer accretion shocks, at ~ 1 — 2i?200, where the run with 
CR feedback produces a stronger SZ signal. The bulk of this effect 
is mostly related to the physical position of the envelope of outer 
accretion shocks in the run with CRs, which is farther from the 
cluster centre of ~ 50 — 100 kpc/h. Again, this is similar to what 
observed in the Zeldovich collapse test (Sect l3.3l l: the physical ef- 
fect of a mixture of CRs and thermal gas in the outer regions pro- 
duces and overall slightly accelerated expansion of the outer shell 
at the same time step. To show that this trend is general to all re- 
simulated clusters, we show in Fig|23]the ratio of the radial pro- 
files of projected maps in clusters HI, H3 and H5, considering both 
X-ray emission (solid lines) and SZ signal (dashed). In order to 
limit the very localized effect of gas clumps with small differences 
in position we post- processed our simulate d maps using a filtering 
technique similar tolRoncarelli et alj |2006), already applied to our 
clusters in lVazza et alj d201 id) . With this technique, the 1 percent 



X-ray brightest pixels at each radial bin are excised from the raw 
maps, and replaced with the average value at the same radius, thus 
highlighting the large-scale behaviour of the smooth gas compo- 
nent only. 

All clusters present a similar deficit of Lx and Y for r < 
0.5i?2O0, of the order of ~ 10 — 30 percent in the case of X-ray 
emission and of ~ 5 — 10 percent in thermal SZ. In the accre- 
tion regions close to ~ -R200, clusters simulated with CR feedback 
tend to have enhanced Lx emission and Y, thus producing slightly 
flatter profiles compared to the runs without CRs . However, a pre- 
cise quantitative estimate is more difficult, due to the systematically 
larger radius of accretion shocks. Taking the most relaxed cluster 
of our sample (HI) as a reference, we find that Lx and SZ effect 
increase at least by ~ 20 — 30 percent. In the case of merging clus- 
ters the maximal difference can be as large as a factor ~ 3 — 4 in 
some projections. 

Combined with the findings of the scaling for (T200, M200) 
and (Lx,20o,A^20o) for the sample of clusters simulated with fixed 
mesh resolution (Sect l4. 1 !2t . these results suggest that the injection 
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and pressure feedback of CRs in large-scale structures can be re- 
sponsible for small but measurable effects. Even if the inclusion 
of additional physical processes in our runs can somewhat alter 
the small-scale distribution of CRs in the innermost cluster regions 
(given that cooling, star formation, AGN and magnetic field should 
be most important for ~ O.lifeoo), our findings for the large scales 
are more robust against the inclusion of these processes, given that 
heating from shock waves is expected to be the most physically 
relevant process at ~ i?200- Interestingly, a recent analysis of a 
large set of clusters observed with ROSAT out to ~ -R200 presents 
the evidence for slightly flatter outer profiles of gas density and 
Lx, compared to the expectations of simple non-radiative simu- 
lated cluster runs dEckert er al, 201 lb). Even if several explanation 
are possible for this, such as gas clumping, the departure from sim- 
ple non-radiative profiles is of the same order as that caused by the 
presence of CRs. 

Hence, the effect of CRs may thus have intriguing con- 
sequences in the physical interpretation of r ecent obse r vation 
of galaxy clusters close to ~ R200 (e.g . iBautz et al 



2009; ISimionescu et alj |201 li lUrban et al 



2009 



2011 



George et al. 

Eckert et al.1 1201 lal) . as well as in the use of cluster scaling re- 
lations to perform "high-precision" cosmological studies (e.g. 
Hall man et alj2006l: IPierre et alj201ll and references therein). 

The application of cosmological simulations including feed- 
back from CRs were performed in the recent past with a modified 
version of the SPH code GA DGET! dSpringell20"05l ; |Pfrommer et al.l 
120071: lJubelgas et all 1 2008). Contrary to our results, these works 
found that the inclusion of CRs leads to an increase of the central 
gas density (and thus the X-ray luminosity and the amplitude of the 
SZ effect) inside ~ O.I-R200, while leaving the outer cluster regions 
essentially un-modified. This is puzzling since the radial profile of 
effecti ve adiabatic index in non-radiative runs of IPfrommer et al] 
d2007l) is very similar to ours (Sect |4.2.T] l. 

Understanding the reason of this difference is not trivial, given 
the different numerical techniques and physical assumptions. We 
argue that most of the differences in the two schemes are due to 
the different way in which the injection of CRs and their transport 
take place. In a recent comparison o f cosmological runs with GAD- 
GET and ENZO dVazza et al]|2011bl) . we investigated in detail how 
matter from accreted gas satellites is shock-heated and distributed 
within the main cluster volume. The result is that while in GADGET 
matter from satellites usually retains its internal low entropy gas, 
delivering it to main cluster centre (R < O.I-R200X m ENZO the 
gas matter from satellites is efficiently shoc k heated and deposited 
at much larger radii, > 0.5i?2oo (see also IVazz3l201 J, for sim- 
ilar results with AM R runs). Very similar c onclusions have been 
recently presented by Simionesc u et al.l d201lh . who compared the 
results on accreted cluster satellites in GADGET and in AREPO. In 
the same work, we also showed that while the total energy flux 
of shock waves in the cosmological volume is quite similar in the 
two codes, the average properties of external accretion shocks are 
not: shocks are sharper, more regular in shape and stronger in grid 
codes. This produces a larger entropy generation in accreted clumps 
in grid simulations, while much more pre-shock entropy genera- 
tion is measured in SPH runs. The effect is much more significant 
than the differences due to simple resolution effects at ~ i?20oin 
the two approaches. Therefore, also the accretion of shock-injected 
CRs from accreted gas clumps would look quite different in the 
two methods, once the injection of CRs is followed in run-time. 
While in GADGET runs we expect that CRs from clumpy accretion 
is delivered in the centre of galaxy clusters, in a grid code such as 



ENZO the injection of CRs (and their dynamical feedback) mostly 
takes place at large cluster radii. The net effect is that while in SPH 
the modification of CRs feedback causes a larger compressibility 
in the centre of structures (and hence an increase of gas density to 
keep pressure equilibrium), in grid codes the modification of CRs 
dominates at large radii, enhancing the compressibility of the outer 
cluster layers already since the early times of structu re formation 
(se e Sectl3.3land Appen dix). Based on our results in Vazza (2011) 
and lVazza et al] d201 lbl). and following previous results from other 



group s (e.g. lAgertz et aU2007l ; IWadslev et ai]|2008l: iMitchell et ail 
|20Q9D, the most likely explanation for the above differences is the 
role of artificial viscosity in the standard SPH formulation, and the 
existence o f sizable pre-shock entropy generation in the standard 
SPH runs dO'Shea etal]|2004l) . Incidentally we notice also that a 
trend for gas density, X-ray luminosit y and gas tempera ture very 
similar to ours has been reported by iRvu & Kand d2004) for grid 
simulations with pressure feedback, where CRs were injected in 
the cluster volume by AGN-like sources. 



5 DISCUSSION 

Our simulations with an implementation of CR physics in ENZO al- 
low us to highlight the role of CRs injected at cosmological shocks. 
Because of the numerical and physical complexity of CRs injection 
and CR feedback on the ICM plasma, a few strong assumptions had 
to be made. 

First, the acceleration efficie ncy at shocks is presen tly uncer- 
tain. We adopted the efficiency o f iKang & Jonesl d2007l) . which is 
based on studies of particles acceleration in SN remnants (how- 
ever in Sec l4.ll we tested also fixed efficiencies). The details of 
particle acceleration in the regime of Mach numbers typical of the 
ICM, M ^ 5 — 10 are not yet robustly constrained by the the- 
ory, due to the difficulty of modelling the large spatial and tem- 
poral scales_jnvolyed_in the diffusive acceleration at such shocks 
(e.g. lKang&Rvull2010h . More recently, interesting studies em- 
ploying particle-in-cells methods investigated additional accelera- 
tion mechanism for particles at shocks (e.g. shock drift accelera- 
tion), su ggesting the possibilities of a different trend with Mach 
number dGargate & SpitkovskvluOllh . Our treatment of CRs injec- 
tion in ENZO runs will enable us to perform cosmological simu- 
lations with any given prescription for the functional dependence 
between acceleration efficiency and background properties of the 
ICM. However, the injection of CRs in the early Universe (z 2), 
and in the most rarefied cosmic environment (p/p CI ,b "C 1) is un- 
certain because the evolution of magneti c fields outside large-scale 
structures is still subject to debate (e.g. IZw eibel & Heitscr] [2008l : 
iMiniati & BelfcOl lHwidrow et alj201 ltlRvu et al]|201 ll and ref- 
erences therein for recent reviews). 

Secondly, we neglect the evolution of the spectrum of acceler- 
ated particles, and their energy losses via Coulomb and/or proton- 
proton collisions. The inclusion of losses may slightly decrease the 
amount of CRs inside our structures. However, for the densities typ- 
ical of the structures of interest here this can be safely neglected, 
given that the time scales for significant energy los ses of CR pro- 
tons are large for typi cal non-cool core clusters (e.g. iMiniatill 20071 : 
IPfrommer etalll2007t) . 

Thirdly, only direct shock acceleration is considered in this 
work as a channel to inject CRs in the cosmic gas. This is a 
safe enough assu mption for most of the sim ulated volume, where 
the turbulent (e.g. iBrunetti & Lazariai feoilbl) and the shock (e.g. 
iBlandford & Eichlenl 19871) re-acceleration of CRs should be neg- 
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Figure 23. Radial average profile of the ratio between X-ray luminosity and SZ effect for re-simulations of clusters HI, H3 and H5 with (V3) and without 
(V0) CR physics. For each cluster, the radial trend for the 3 projections are shown in different colours. 



ligible for the usually simple geometry of large-scale accretion 
pattern. In the centre of clusters, however, turbulence and shock- 
lets structures can introduce also some significant re-acceleration 
of CRs. This additional channel of acceleration is not expected 
to boost CR energy ratio by more than a factor ~ a few over ~ 
Gyr time-scales. However, a future development of our method 
will also include these effects (and their feedback on the gas en- 
ergy) in run-time. The injection of CRs at supernovae and AGN is 
also neglected in our model. Numerical studies using GADGET sug- 
gested that these additional CR sources are relev ant only on very 
small ~ 10 kpc scales in proximity of the sources ( Pfrommer et al. 
2007). Heating of the therm al gas from the Alfven waves excited 
by streaming of CRs (e.g. IGuo & Ohll2008l) is also an interesting 
physical process which can be incorporated in future work. 

The description of the ICM in the runs presented here is sim- 
plified, since we neglected radiative cooling, and feedback from 
star formation and galactic activities, and magnetic fields. Radia- 
tive cooling can have an impact on t he injection of CR en ergy at 
cold lumps of gas in the ICM (e.g. IPfrommer etafl 12007). How- 
ever, the handling of cooling in simulations is still problematic 
since it produces strong (and un-observed) cooling flows in the in- 
ner regions of clusters, and it is still unclear how to quench them 
self-consistently. Magnetic field, even if not dynamical important, 
can be relevant to the study of CR energ y in clusters because it 
affects the spatial pr opagation of CR (e.g. lHanasz & Lescbll2000l : 
lJubelgas etafl l2008). However, for typical conditions in the ICM 
the scale-length of CR diffusion is smalle r, or very close to , than our 
highest spatial resolution, 25 kpc/h (e.g. lBlasi et al ] |200l and ref- 
erences therein), making our results still valid in a statistical sense. 

Within the range of assumptions listed above, our method pro- 
duces a number of interesting results with consequences on the 
modelling of galaxy clusters, and their comparison with observa- 
tional data. Our runs at fixed mesh resolution show that the largest 
effect of CR feedback is concentrated at the outer accretion re- 
gions of cosmic structures, at over-densities close to the critical 
one (Sect |4. it . At p/p cr ,b ~ 1 we measure P cr /P s ~ 0.3 — 0.5 
on average. This result is very stable against all explored modifi- 
cations of numerical parameters (the pressure ratio between CRs 
and gas for p/p CI ,b ^ 1 is, on the other hand, significantly de- 
pendent on the assumed parameters). The use of AMR enables the 
investigation of the innermost region of galaxy cluster at a higher 
resolution (25 kpc/h). In this work we analysed 7 galaxy clus- 
ters with different dynamical states, in the range of total masses 



2 ■ 1O 14 M < M < 1.1 ■ 10 1S M© (Sectl^a. We report a very 
low pressure ratio in the innermost cluster regions (^ O.I-R200) 
is very small, P CI /P g ~ 5 ■ 10~ 3 — 3 • 10~ 2 , increasing up to a 
~ 0.1 at i?200- No obvious dependence on cluster mass or dynam- 
ical state can be found in this small sample, except for the trend 
that dynamically perturbed systems display slightly larger values of 
Pcr/Pg, and tend to concentrate a significant amount of CR pres- 
sure in clumps of accreted sub-halos. In all systems, the effective 
adiabatic index of the mixture of gas and CRs is very close to the 
mono-atomic value of 7 = 5/3 for ^ 0.5i?2oo- In a few perturbed 
systems, however, it can be as small as y e g ~ 1.6 close to i?20o- 

It is important to compare the amount of CRs in our simulated 
clusters with available data from observations. A direct approach to 
constrain the energy content of CR protons in galaxy clusters is the 
observation of 7-ray emission from the decay of the neutral pions 
due to proton-proton collisions in the ICM. Gamma ray upper lim- 
its from EGRET observations allo w to put limits of E CY /E K < 0.3 
in several nearby galaxy clusters dReimer et al]|2003l) . More strin- 
gent limits are derived from deep pointed observations at energies 
of > 100 GeV with Cherenkov telescopes. These limits depend 
on the (unknown) spectral shape of the proton-energy distribution. 
For the relevant case 5 = 2 (with N(p) oc p~ s ) the limi ts are 
E CI /E g < 0.1 dAharonian etai]|200a [Xleksic et alj|201fj|) . and 
they are less stringent for steeper spectra. Recently, the FERMI 
satellite greatly improved the sensitivity of observations at GeV 
energies, and the present upper limits for a large sample of nearby 
galaxy clusters are _E cr /_E B < 0.05 (e.g. [Ackerman n"et al.ll201u 
I Jeltema & Profumol I201TI) . with a poor dependence on S. In ad- 
dition to these methods, also the limits to the presence of diffuse 
Mpc-scale radio emission in clusters can be used to constrain sec- 
ondary electrons and thus the energy density o f the primary CR 
protons terunetti et"aT]|2007l ; iBrown et al]|201ll) . In this case, the 
limits depend also on the cluster magnetic field strength and are 
complementary to those obtained from 7-ray s. In the relevant case 
of an average magnetic field in cluster of a ~ pG, radio obser- 
vations of clusters with no Mpc-scale radio emission suggest that 
E CI /E S ^ 0.05, while the limits are less stringent for smaller aver- 
age magnetic fields. These limits usually refer to innermost ~ Mpc 
regions of clusters, where both the number density of thermal pro- 
tons and the magnetic field are larger. At present no tight constraints 
are available for the clusters outskirts, where the CR contribution 
might be larger. 

Given the mean value measured inside the innermost regions 
of all our re-simulated clusters (Sect l4.2l (. we conclude that at 
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present no obvious tension exists between the estimated energy 
budget of CRs injected at cosmological shock wav es, at least for 
accele ration efficiencies in the range of those of Kang & Jonesl 
d2007l) . It is still possible that the inclusion of additional sources 
of CRs (as supernovae and AGN), as well as of re-acceleration 
from turbulence and shocks may increase the pressure of CRs close 
to available upper-limits. Therefore, in the next future a detailed 
comparison of 7-ray fluxes from secondary particles (possibly with 
more CR acceleration mechanisms) and observations will be im- 
portant to address this topic more robustly. 

When we compare the 3-D structure of the thermal gas in 
these clusters to their re-simulations which only employ standard 
non-radiative physics (Sec l4.2. It . we find small systematic trends: 
the innermost gas density, gas pressure and entropy at z=0 are a 
few percent smaller if feedback from CRs is considered, while they 
are larger by a similar amount at ~ R200, leading to slightly flat- 
ter radial profiles. Based on our 1-D tests with the Zeldovich col- 
lapse (Sect |3.3t and also on specific re-simulations discussed in the 
Appendix A, we conclude that this arises from the time-integrated 
effect of having a lower 7 e ff at the outer accretion regions of clus- 
ters. Accretion shocks are strong, M ^> 10, and their acceleration 
efficiency is large. CRs are accumulated here and they become dy- 
namically important as the outer shells of matter of the structures 
expands with the growth of structures, because of their softer equa- 
tion of state. This enhances the total pressure jump felt by infalling 
gas at ~ -R200 (while the gas pressure jump is slightly decreased 
by the effect of having a reduced thermalization efficiency), it in- 
creases the effective Mach number of the outer shocks, and it causes 
a slightly faster expansion of the outer shells compared to standard 
non-radiative runs (Sec |3.3t . As a net effect, a few percent more gas 
matter is deposited in the outer regions, and a slightly less dense 
core is formed inside clusters. The reduced thermalization and the 
decreased gravitational potential in the innermost cluster regions 
also cause a smaller gas pressure, compared to simulations without 
CR feedback, by a few percent. 

As we discussed in Sect |4.2.2l this feature has significant ef- 
fects on the simulated properties of X-ray emission and SZ signal 
for our sample of clusters. When re-simulations with and without 
feedback from CR are compared, both signals are found to be de- 
creased by a 10 — 30 percent in the innermost cluster regions (^ 
O.I-R200) if CRs are accelerated, while they are increased by factors 
~ 0.5 — 4 at i?200- The reported trends are found to be at varianc e 
with other ones based on SPH simulations (Pfr ommer et alj|2007h . 
where a significant boost of X-ray and SZ signal is reported for the 
innermost cluster regions. To our understanding these differences 
can be totally ascribed to the fundamental differences in the way 
the two numerical approaches model at a run-time the accretion of 
gas satellites, their entropy enrichment and their stripping into the 
main cluster atmosphere. This problem is related to the more funda- 
mental ones of artificial viscosi ty and pre-shock entropy generation 
in SPH (e.g.lAgertz et alj2007l;IWadslev et al.l2008l:lMitchell et all 
l2009l;IVazzall201 ll i lVazza et alfeoi lbl i lSiiacki et ai]|201l[ and ref- 
erences therein). 



6 CONCLUSIONS 

In this work we presented our numerical implementation of cos- 
mic ray injection, advection and feedback in cosmological simula- 
tions with the ENZO code. Our method incorporates a prescription 
for Diffusive Shock Acceleration (e.g. lBelllll978l ; iDrurv & Voelkl 



ll98ll ; lKang& Jonesll2007t) . and channels at run-time energy from 
the thermal gas to the CR pool, thus reducing the post-shock ther- 
mal energy flux. The main step of our algorithm are the following: 
a) in run-time we measure the 3-D distribution of Mach number us- 
ing a shock finder based on pressure jumps; b) we estimate the total 
energy flux of CR protons as a function of Mach number, assum- 
ing an acceleration e fficiency from theoretical models (in our case, 
Kan g & Jonell2007l) : c) we update the energy of CR energy within 
each shocked cell, and we advect CR energy assuming no diffusion 
of CRs, for a fixed adiabatic index of 7 cr = 4/3; d) we feed the 
total pressure, P g + P cr (rather than the simple gas pressure, as in 
standard runs) into the Riemann solver. Therefore, the composite 
fluid of gas and CRs in the simulation follows an effective adia- 
batic index, 7 e ft ^ 7 = 5/3, depending on the local energy ratio 
between CRs and the gas energy within the cells. 

We tested our method in simple 1-D tests (Sec(3](, finding 
good agreement with analytical estimates for shock-tube tests. The 
role of the various parameters involved in the injection and advec- 
tion of CRs was tested separately with shocks of different strength 
(Sec l3.2l >. The cosmological application of our method for CRs was 
studied in 1-D using the Zeldovich pancake collapse test, where 
also the use of adaptive mesh refinement was tested (Sec l3.3b . 

We studied the injection and the evolution of CRs in large- 
scale structures with cosmological simulations at fixed grid reso- 
lution CSec l4. lb and with adaptive mesh refinement (Sec l4.2t . For 
fixed grid runs, we report an increase in the total CPU time of the 
order of ~ 10 — 20 percent compared to standard non-radiative 
runs. These low and moderate resolution runs were used to esti- 
mate the large-scale properties of CRs in the cosmic volume, and 
to test the robustness of our findings against a number of assump- 
tions. The level of CR energy inside cosmic structures is found to 
be small, P cr /P g ^ 0.1, with a peak at the over-density typical of 
outer accretion regions. We report that only the distribution of CRs 
outside of cosmic structures is strongly dependent on the details 
involved in the acceleration in the early cosmic epochs (and in the 
most rarefied environments) while the distributions of CRs are very 
stable for the innermost regions of clusters. 

Using adaptive mesh refinement we investigated the properties 
of 7 galaxy clusters, with masses in the range 2.1 • 1O 14 M0 ^ 
M ^ 1.5 ■ 10 15 A/0, and different dynamical states. The dynamical 
role of accelerated CR energy is always quite small, and plays a 
significant dynamical role only close to ~ i?200- In the centre of 
clusters instead the pressure of CRs is small, P cr /P g ~ 0.02 — 
0.05. These values are presently consistent with the upper-limits 
provided from 7-ray observations (Sec[5]l. 

The effects of CRs on the overall evolution of clusters have 
small and systematic effects on the 3-D distribution of the ther- 
mal baryonic gas. In all re-simulated clusters in the innermost 
gas density, temperature and entropy are reduced by a few per- 
cent, while they are enhanced on average by the same amount at 
P200 (Sec 14.21 . This comes from the fact that CRs first modify the 
compressibility of outer accretion regions during the formation of 
structures, leading to an enhanced post-shock compression and to 
a slightly faster expansion of the outer cluster layers compared to 
standard simulations. This produces also a corresponding decrease 
of X-ray emission and of the thermal SZ signal from the inner clus- 
ter region, and an enhancement of a factor ~ 0.4 — 4 close to R200, 
and depending of the dynamical state of the clusters (Sec 14.2.21 . 

These systematic trends in galaxy clusters are at variance with 
results obt ained with SPH for qua l itatively similar treat ment of CRs 
dynamics (Pfro mmer et alj|2007l ; ljubelgas et alj|2008l) . where en- 
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hanced gas density is found in cluster cores, and unchanged gas 
distributions are measured further out. The reason of these differ- 
ences is not clear. We conjecture that the most likely explanation 
lies in the different way in which transport and mixing motions are 
modelled in grid methods and in SPH, leading to a more efficient 
large-scale mixing of injected CRs in Eulerian runs. This would be 
in line with recent comparisons performed by our group on simple 
non-r adiative cosmolog ical simulations using and ENZO and GAD- 
GET dVazza et"al]|201 lbl) . Performing similar comparisons in cos- 
mology, and with the inclusion of CRs physics, is a necessary next 
step to study non-thermal processes in cosmology. 

To summarize, the first step for incorporating particle accel- 
eration at shocks in run-time in ENZO, and the physical feedback 
of CRs on the baryon gas is successful for the ensemble of tests 
and cosmological simulations we explored. To our knowledge, this 
is the first time that such studies are successfully performed for 
cosmological simulations with adaptive mesh refinement, and for 
realistic models of shock acceleration and reduced thermalization. 
This model enables us to explore non-thermal phenomena in galaxy 
clusters, and to take advantage of the capabilities of the PPM 
method to model turbulence and shocks in the simulated intra clus- 
ter medium. The inclusion of more co mplex treatmen t of CR dy- 
namics (e.g. using spectral energy bins, Miniati 2007), as well as 
the porting of these meth ods to more sophist icated and recent de- 
velopment of ENZO (e.g. ICollins et alj|201oh are under way, and 
will allow us to simulate the non-thermal features of the ICM with 
unprecedented detail. 
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APPENDIX A: INVESTIGATING THE REASON OF A 
LOWER CENTRAL DENSITY 

In this Section we present a numerical test designed to highlight 
the impact of CR physics on the formation of galaxy clusters, at 
an early stage of their assembly. In this test, we selected a cu- 
bic sub- volume of side 10 Mpc/h around the forming cluster H3 
at z=2. This is a standard non-radiative run in the redshift range 
30 ^ z > 2. Starting from z — 2, we evolve two re-simulations 
of the same volume, allowing CRs injection and feedback in the 
first, and following standard non-radiative physics in the second. 
At each root-grid time step (At ~ 0.05 Myr) we compare the 
thermal properties of gas halos in the two runs. This highlights the 
timely role of CRs in our runs, which is responsible of the signif- 
icant differences discussed in Sect l4.2l In the first three panels of 
Fig lAll we show the time-sequence of profiles for the normalized 
difference of gas density, gas pressure and total pressure for a mas- 
sive halo in our sample (M ~ 4 ■ 1O 13 M0 at z=2). In detail, we 
compute (Vcr{R) - Vo(R))/Vo(R), where V CI (R) is the profile of 
each quantity in the run with CRs feedback, and Vo(R) is the cor- 
responding profile in the run with no CRs feedback. The last three 
panels of the same Figure report the time evolution of CR pressure, 
pressure ratio and 7 c g for the re-simulation with CRs feedback. 

The injection of CRs within halos starts mainly from outer 
accretion shocks (at ~ 0.2 Mpc/h in the Figure), and at internal 
merger shocks at the centre of the structure. The energy feedback 
of CRs is initially more important in the outer proto-cluster regions 
(since the Mach number is large here), leading to a small and con- 
tinuous decrease of the effective adiabatic index of the mixture of 
gas and CRs. The gas pressure in the run with CRs after a few 
time-steps is initially reduced by a ~ 1 — 2 percent, while the total 
pressure remains almost constant. The post-shock gas density at the 
outer regions is progressively increased by the enhanced compress- 
ibility at the shock. This in turn leads to a smaller deposition of gas 
matter (reducing gas density by ~ 1 percent) at the outer cluster 
parts, compared to the standard run. The outer shells expand into 
the rarefying cosmic volume, and the smaller adiabatic index of CR 
energy makes it even more dynamically important as the expansion 
goes on. and By the end of our test (z ~ 1.6) this has caused the 
increase of ~ 5 percent in the total (gas+CR) pressure at the outer 
proto-cluster regions. We thus expect that shock waves are slightly 
stronger at ~ R200 in the run with CRs, even at later times. This 
makes the injection of new CRs an even more efficient process, and 
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Figure Al. First three panels: time evolution of the ratio of gas density, gas pressure and total pressure in a re-simulated halos of run H3 with and without CR 
physics. Last three panels: time evolution of the CRs pressure, pressure ratio and effective adiabatic index for the same cluster. The injection of CR has been 
started at z=2 in this case. The time sequence goes from blue to green to red colours, the time lag between each profile is At ~ 0.05 Gyr. 



maintains the difference with standard runs even at later times, as 
observed in our clusters at z=0 (Sect |4.2| l. 



